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Abstract 

Very little is known about genetic factors that regulate life history transitions during ontogeny. Closely related tiger salamanders 
(Ambystoma species complex) show extreme variation in metamorphic timing, with some species foregoing metamorphosis alto- 
gether, an adaptive trait called paedomorphosis. Previous studies identified a major effect quantitative trait locus (metl) for meta- 
morphic timing and expression of paedomorphosis in hybrid crosses between the biphasic Eastern tiger salamander (Ambystoma 
tigrinum tighnum) and the paedomorphic Mexican axolotl (Ambystoma mexicanum). We used existing hybrid mapping panels and a 
newly created hybrid cross to map the met 7 genomic region and determine the effect of metl on larval growth, metamorphic timing, 
and gene expression in the brain. We show that met 7 maps to the position of a urodele-specific chromosome rearrangement on 
linkage group 2 that uniquely brought functionally associated genes into linkage. Furthermore, we found that more than 200 genes 
were differentially expressed during larval development as a function of met 7 genotype. This list of differentially expressed genes is 
enriched for proteins that function in the mitochondria, providing evidence of a link between metl , thyroid hormone signaling, and 
mitochondrial energetics associated with metamorphosis. Finally, we found that metl significantly affected metamorphic timing in 
hybrids, but not early larval growth rate. Collectively, our results show that metl regulates species and morph-specific patterns of 
brain transcription and life history variation. 
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Introduction 

Some closely related species exhibit dramatic phenotypic dif- 
ferences of adaptive significance. Such species provide excel- 
lent models to investigate the genetic basis of complex traits 
and life history evolution. An important and difficult objective 
in working with natural species is to identify gene and gene 
functions that are associated with phenotypic differences 
(Jones 1998; Bradshaw and Schemske 2003; Hoekstra et al. 
2006). In some cases, potential candidate genes are predict- 
able given the nature of the phenotypic differences, such as in 
the case of color polymorphisms (Mundy 2007) and craniofa- 
cial variations in birds and fishes (Schneider and Helms 2003; 
Streelman et al. 2007). However, for most traits, it is difficult 
to predict genetic and developmental components of the se- 
lection process. When it is possible to cross species and seg- 
regate molecular and phenotypic variation, quantitative trait 



locus (QTL) mapping provides an unbiased and efficient 
method to initially locate genetic factors in the genome. 
Subsequent comparative and functional genomic approaches 
can then be used to narrow the search for causative genetic 
factors and reveal insights about gene functions, develop- 
ment, and evolution. 

Here, we consider genetic factors that activate postembry- 
onic developmental programs during ontogeny to enable crit- 
ical life history transitions. Amphibian and insect larvae 
undergo a dramatic life history transition known as metamor- 
phosis, during which adult traits develop to enable strategies 
for predator avoidance, dispersal, resource utilization, and 
mating (Wilbur 1980). Thyroid hormone (TH) is the master 
regulator of amphibian metamorphosis, and release of TH 
by the thyroid is under complex control of the hypotha- 
lamic— pituitary— thyroid axis (Shi 2000; Denver et al. 2002). 
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At a molecular level, the effects of TH are mediated by TH 
receptors, which are nuclear transcription factors that regulate 
the transcription of genes in a TH-dependent manner 
(Buchholz et al. 2006). Thus, rising TH levels during larval de- 
velopment bring about metamorphosis by triggering tran- 
scriptional programs that activate tissue remodeling, cellular 
differentiation, and tissue regression — developmental pro- 
cesses that transform larval phenotypes into adult phenotypes 
(Shi 2000). According to this model, variation in metamorphic 
timing may trace to variation in the timing of TH-dependent 
transcriptional events during development. 

The timing of metamorphosis is a critical life history trait 
that evolves in response to aquatic habitat characteristics that 
determine the rate and plasticity of larval growth and the 
length of the larval period (Wilbur and Collins 1973; Wilbur 
1980; Semlitsch and Wilbur 1989). Species that use relatively 
permanent habitats for larval development generally have 
longer larval periods than species that use ephemeral habitats. 
For example, some anuran tadpoles hatch in rapidly drying 
desert ponds and metamorphose after only 8 days of devel- 
opment (Newman 1989), whereas salamanders in perma- 
nent, high altitude lakes may overwinter as aquatic larvae 
and transform the following year, or forgo metamorphosis 
altogether (Petranka 1 998). The observation of life history var- 
iation between anuran and salamander species is not surpris- 
ing because these groups diverged from a common ancestor 
more than 300 Ma. However, there are also substantial dif- 
ferences in developmental timing within both of these orders, 
of which the pronounced differences in metamorphic timing 
between closely related tiger salamander lineages 
(Ambystoma spp.) are an unusually dramatic example. Such 
species groups provide models for identifying genes for meta- 
morphic timing and investigating the relationships between 
development, life history variation, and life cycle evolution 
(Voss et al. 2003, 2012; Voss and Smith 2005). 

In most Eastern tiger salamander (A. tighnum tigrinum) 
populations, larvae undergo an obligate metamorphosis 
after a variable period of growth (Petranka 1998). In contrast, 
Mexican axolotls (A. mexicanum) do not produce a high titer 
of TH during larval development and thus do not metamor- 
phose (Prahlad 1968; Galton 1992), a process called paedo- 
morphosis (Gould 1 977). As a consequence, sexually mature 
A. mexicanum retain juvenile traits and complete their life 
cycles in the aquatic habitat. Crosses between A. t. tighnum 
and A. mexicanum, followed by backcrosses of F1 hybrids to 
A. mexicanum, show that a major effect QTL (metl) explains 
both continuous variation in metamorphic timing and expres- 
sion of paedomorphosis (Voss and Shaffer 1997, 2000; Voss 
and Smith 2005; Voss et al. 2012). Heterozygotes that inherit 
a single A. t. tigrinum metl allele almost invariably undergo 
metamorphosis. However, homozygotes that inherit two 
A. mexicanum metl alleles metamorphose at a later time 
than heterozygotes, or fail to undergo metamorphosis alto- 
gether and exhibit paedomorphosis. The number of 



paedomorphs observed from hybrid crosses depends on ge- 
netic background. For example, when hybrid crosses are made 
using A. mexicanum from a laboratory population, the ob- 
served number of paedomorphs is statistically consistent 
with a single gene model (Voss 1995). However, when 
hybrid crosses are made using wild-caught A. mexicanum, a 
more complicated genetic basis is implicated and fewer pae- 
domorphs are observed than expected under a single gene 
model. In these crosses, most of the A. mexicanum metl ho- 
mozygotes do not exhibit paedomorphosis. Instead, they 
metamorphose on average 36 days later than siblings carrying 
an A. t. tigrinum metl allele (Voss and Shaffer 1997, 2000; 
Voss and Smith 2005). Thus, met 1 alleles from paedomorphic 
A. mexicanum delay metamorphosis while metl alleles from 
metamorphic A. t. tigrinum decrease the time to metamor- 
phosis, and this consequently decreases the probability of 
expressing paedomorphosis. 

In this study, we used comparative, quantitative, and func- 
tional genomic approaches to map the metl genomic region 
and investigate the effect of met 1 genotype on larval growth, 
neurodevelopmental gene expression, and metamorphic 
timing. Our previous work identified metl using anonymous 
molecular markers (Voss and Shaffer 1997, 2000). Subse- 
quently, an expressed sequence tag (EST) (ctg325) was iden- 
tified for metl that shows nucleotide similarity to human ngfr 
(Voss and Smith 2005). We show here that ngfr and metl 
map to a genomic region that was structured by chromosomal 
rearrangements after the divergence of salamanders and an- 
urans from a common ancestor. As a result, a unique set of 
metl region genes were brought together during evolution, 
and several of these encode proteins that function in related 
biological processes. Further, several genes in linkage disequi- 
librium with the metl region are differentially expressed in the 
brain during larval development, indicating that more than 
one gene from the metl region may explain phenotypic dif- 
ferences between these species. At a genomic level, more 
than 200 genes, including genes that code for proteins with 
mitochondrial functions, were differentially expressed be- 
tween metl genotypes that significantly affect metamorphic 
timing. These changes were observed during a specific interval 
of time late in the larval period in interspecific hybrids, but at a 
much earlier time in the parental species. Our results show 
that metl regulates temporal changes in brain transcription 
during larval development and the timing of metamorphosis, 
but not the rate of larval growth preceding these events. 

Materials and Methods 

Salamander Care and Genotyping 

Second-generation hybrids were generated by backcrossing 
an A. t. tigrinum-A. mexicanum Ft hybrid male to an A. mex- 
icanum female via in vitro fertilization (Voss 1995). Embryos 
were transferred to large containers of pond water and 
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aerated. Hatching was synchronized manually and from 
hatching onward, larvae were reared individually. From hatch- 
ing until approximately 21 days posthatching (DPH), larvae 
were fed brine shrimp napuli (Artemia). After this, larvae 
were fed California blackworms (Lumbriculus). At 19 DPH, 
tail clippings were collected for genomic DNA extraction 
using phenol/chloroform. Genotyping was performed to 
determine whether individuals were carrying A. mexicanum, 
A. t. tigrinum, or a combination of species-specific alleles at 
loci from the metl genomic region. Genotypes for several loci 
(discussed later) were obtained using AcycloPrime-FP chemis- 
try and a Wallac Victor3 plate reader as described in Smith 
et al. (2005). The polymerase chain reaction primers and ex- 
tension probes that were used to type metl region loci are 
presented in supplementary table S1, Supplementary Material 
online. 

At 28, 42, 56, 70, 84, 98, and 1 12 DPH, 18 unique indi- 
viduals were measured for snout-vent length (SVL) prior to 
tissue harvesting (i.e., no individual was measured more 
than once). Whole brains and pituitaries were then removed 
and immediately flash frozen in liquid nitrogen (six of each 
metl genotype per time point). Excluding two recombinants 
(one at 70 DPH and one at 1 1 2 DPH), half of these individuals 
were homozygous for axolotl alleles across ngfr, rasdl, 
map2k3, rail, and shmtl (i.e., had a m et1 mexlmex genotype), 
whereas the other half were heterozygous for these loci (i.e., 
had a metl mex/Att genotype). An additional 30 animals, includ- 
ing one recombinant {jygff™"**, rasdl mex/mex , map2k3 mex/Att , 
rail"* 3 *"*, and shmtl mex/mex ), were allowed to develop indef- 
initely, so that we could examine the relationship between 
metl genotype and metamorphic timing. All animal care 
and use procedures were in accordance with protocols ap- 
proved by the University of Kentucky Institutional Animal 
Care and Use Committee. 

Comparative and QTL Mapping of metl 

Previously, we performed A. mexicanum/A. t. tigrinum x 
A. mexicanum backcrosses and developed markers from 
ESTs to make linkage maps (Voss 1995; Voss and Shaffer 
1997, 2000; Voss and Smith 2005; Voss et al. 201 1). In this 
study, we used the mapping panel and markers from Voss 
et al. (201 1) to perform genome-wide QTL scans for meta- 
morphic timing QTL. Genome-wide marker information is de- 
tailed in Voss et al. (2011). Marker orders were determined 
using Multipoint 2.2 (MultiQTL Ltd., Hafia, Israel) and the 
Kosambi mapping function (Kosambi 1944). QTL were iden- 
tified using R/qtl (Broman et al. 2003), the scanone function, 
the normal QTL model, and Haley-Knott regression (Haley and 
Knot 1992). Genome-wide thresholds for evaluating the sig- 
nificance of QTL were determined from 1 ,000 replicated data 
sets (Churchill and Doerge 1994). We also used the mapping 
panel from Voss and Smith (2005) to map additional loci and 
more finely map the position of metl. Primer sequences, 



diagnostic polymorphisms, and polymorphism detection 
assays for metl region loci are summarized in supplementary 
table S1, Supplementary Material online. 

Statistical Analyses of Growth and Metamorphic Timing 
in Backcrossed Hybrids 

We assessed whether backcrossed hybrids with nonrecombi- 
nant heterozygous and homozygous genotypes for ngfr, 
rasdl, map2k3, rail , and shmtl differed in growth trajectory 
by fitting a general linear model of the following form: 
SVL t ^= p t + G t + T i+ (GT) e -+ T 2 i+ (GT 2 ) 

ti+£tij, where p t is the 
intercept for metl mex/mex individuals, G t is the difference be- 
tween the intercept for met1 mex/mex and metl mex/Att individ- 
uals, Tj is the linear regression coefficient for metl mex/mex 
individuals, (GT) t/ is the difference between the linear regres- 
sion coefficients for me tl mexlmex and metl mex/Att individuals, 
T 2 ; is the quadratic regression coefficient for metl mex/mex indi- 
viduals, (GT 2 ) ft is the difference between the quadratic regres- 
sion coefficients for metl mex/mex and metl mex/Att individuals, 
and E ti j is the error term for the y'th individual of genotype f 
(a dummy variable) sampled at time /'. This error term is as- 
sumed to be normally distributed with mean = 0 and vari- 
ance = a 2 and describes the deviation of each individual 
(irrespective of genotype) from its expected value based on 
the fitted model, which accounts for systemic differences be- 
tween the genotypes when making predictions. A backward 
selection procedure was then applied in which parameters 
with (--statistic P values > 0.05 and deletion test P values > 
0.05 were removed from the model (Crawley 2007). When 
conducting this procedure, the quadratic interaction term be- 
tween genotype and time (i.e., (GT 2 ) t/ ) was removed prior to 
the linear interaction term between genotype and time (i.e., 
(GT) fi ), which was removed prior to the main effect of geno- 
type (i.e., G t ). We also assessed whether biphasic salamanders 
that were nonrecombinant for ngfr, rasdl , map2k3, rail , and 
shmtl differed in time to completion of metamorphosis by 
conducting a one-tailed Mann-Whitney test. Salamanders 
that did not complete metamorphosis within 365 DPH were 
scored as paedomorphs. 

RNA Extraction, Microarray Platform, and Quality Control 

Total RNA was extracted individually from the brains of back- 
crossed hybrids sampled between 42 and 112 DPH. RNA ex- 
tractions were conducted using TRIzol (Invitrogen) according 
to the manufacturer's protocol. Samples were further purified 
using RNeasy mini columns (Qiagen). Upon extraction, RNA 
samples were quantified via UV spectrophotometry and qual- 
ified via an Agilent Bioanalyzer. Four RNA samples per geno- 
type (i.e., nonrecombinant homozygotes vs. heterozygotes) 
from each of six time points (42, 56, 72, 84, 98, and 112 
DPH) were submitted to the University of Kentucky 
Microarray Core Facility (UKMCF). The UKMCF generated 
biotin-labeled cRNA targets for each of the 48 samples and 
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independently hybridized each sample to a custom 
Ambystoma Affymetrix GeneChip (Page et al. 2007). This 
microarray platform consists of 4,844 probe sets that were 
designed using a preferred maximum of 1 1 probe pairs 
(match/mismatch) and a possible minimum of 8 probe pairs. 
Of these 4,844 probe sets, 57 are Affymetrix controls, 4,596 
were designed from A. mexicanum sequence, and 191 were 
designed from A. t. tigrinum sequence. We subjected all mi- 
croarray data to low-level quality control (QC) as described 
previously (Page et al. 2010) and used the robust multiarray 
average (RMA) algorithm (Irizarry et al. 2003) to generate 
probe set-level expression summaries. We then examined 
the probe set-level data by computing correlation matrices 
for replicate GeneChips, conducting principal component 
analysis, and rendering pair-wise M versus A plots of replicate 
GeneChips. Upon conducting low-level and probe set- 
level QC, it was clear that one of the met i mex/mex 42 
DPH GeneChips was aberrant. Thus, this GeneChip was 
removed prior to further processing and analysis. All QC 
analyses were performed using the R statistical comput- 
ing environment (www.r-project.org, last accessed July 16, 
2013) and the affy and affyPLM packages (Gautier et al. 
2004; Bolstad et al. 2005). All sequence data, the raw 
microarray data (.CEL files), and the microarray annotations 
are available at Sal-Site (http://www.ambystoma.org/genome- 
resources, last accessed September 17, 2013). 

Identification of Zero Mismatch Probe Sets and 
Expression Summarization 

Despite the shallow phylogenetic distance between A. mex- 
icanum and A. t. tigrinum (Shaffer and McKnight 1996; 
O'Neill et al. 201 3), it is possible that targets from backcrossed 
hybrids will have variable hybridization efficiencies for probes 
on the Ambystoma GeneChip due to sequence divergence 
between the parent species (Bar-Or et al. 2007; Buckley 
2007). To address this issue, we used probe sequences as 
queries in Blast searches of A. mexicanum and A. t. tigrinum 
EST contig databases. Blast alignments were then used to ex- 
trapolate the number of mismatches between each probe's 
sequence and A. mexicanum and A. t. tigrinum EST-based 
contigs, respectively (Page et al. 2010). In total, 1,604 probe 
sets (-33%) were designed from contigs with predicted 
orthologs in A. mexicanum and A. t. tigrinum, and for 1 ,320 
(-82%) of these, we did not detect any mismatches between 
the evaluated probe sequences and the contig sequence 
from the heterologous species. However, for 31 of these 
probe sets (-2%), we were only able to assess a subset of 
the probes, and we therefore refer to these probe sets 
as equivocal zero mismatch probe sets (EZMMPSs) to distin- 
guish them from the 1,289 zero mismatch probe sets 
(ZMMPSs) that were fully evaluated. Upon completing this 
exercise, we used the subset option of the RMA function 
in the affy package to calculate probe set-level expression 



summaries based on the 1,320 ZMMPSs and EZMMPSs and 
used this RMA matrix to identify differentially expressed genes 
(DEGs; discussed later). 

Statistical Analysis of Microarray Data 

To arrive at a list of DEGs, we implemented two statistical 
approaches. First, the limma package (Smyth 2004) was 
used to compare the alternate metl genotypes at 42, 56, 
70, 84, 98, and 112 DPH. This involved fitting a linear 
model to each probe set that contained 12 coefficients, one 
for each time by genotype combination. We then reoriented 
this model in terms of contrasts between the two genotypes 
at each of the six time points. We corrected for multiple test- 
ing by applying a 0.05 false discovery rate (Benjamini and 
Hochberg 1 995) to the P values associated with the moder- 
ated F statistics that resulted from our contrast matrix. In 
addition, for some downstream analyses, we required that 
DEGs identified using limma be differentially regulated by 
> 1 .5-fold as a function of metl genotype at one or more 
time points. 

To identify genes with modest but temporally consistent 
expression differences between the metl genotypes, we also 
used the maSigPro package to conduct a stepwise qua- 
dratic regression analysis (Conesa et al. 2006). We corrected 
for multiple testing by applying a 0.05 false discovery rate 
(Benjamini and Hochberg 1995) at the overall model level. 
We then implemented a backward selection procedure 
(Conesa et al. 2006) that removed parameters from the full 
quadratic model with P values > 0.001, a stringent threshold 
that better enabled us to exclude genes with modestly differ- 
ent temporal trajectories between metl genotypes, but lit- 
tle to no difference in elevation (i.e., expression levels). We 
also required that models have /? 2 >0.50 and significant ge- 
notype, genotype by time, or genotype by time 2 terms. Finally, 
we graphically inspected the 96 genes that met the above 
criteria and manually selected genes that were robustly differ- 
ent between the alternate metl genotypes at several time 
points. 

Enrichment Analyses 

We used EASE analyses, as implemented by DAVID (Dennis 
et al. 2003), to identify biological process, cellular component, 
and molecular function gene ontology terms that were statis- 
tically enriched in lists of DEGs. This approach applies a con- 
servative adjustment to the calculation of Fisher's exact 
probability from a hypergeometric distribution (Hosack et al. 
2003). The ZMMPSs and EZMMPs with predicted orthologies 
to human proteins were used to generate expected values 
(i.e., as the background). The count threshold was set 
to two, and the EASE threshold (i.e., the critical value of the 
corrected Fisher's exact P value) was set to 0.05. 
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Examination of Hybrid DEG Expression in the Parent 
Species 

In a previous study, we used the same microarray platform 
used here to investigate neurodevelopmental gene expression 
in A. mexicanum and A. t. tigrinum at 42, 56, 70, and 84 DPH 
(Page et al. 2010). We used the raw data from this prior ex- 
periment (i.e., the .CEL files) to generate a zero mismatch 
RMA matrix for the parent species that is analogous to the 
zero mismatch RMA matrix for backcrossed hybrids described 
earlier. We then obtained log 2 (RMA) values from this matrix 
for probe sets that were identified as differentially expressed 



between mef 7 



and metf 



hybrids in this study 



using limma. For each of these probe sets, we used 
Welch's f statistic to test for differential expression between 
A. mexicanum and A. t. tigrinum at each time point that we 
sampled parental gene expression (Page et al. 2010). In addi- 
tion, we used Kruskal's nonmetric multidimensional scaling 
(NMDS) to examine relationships among hybrid, A. mexica- 
num, and A. t. tigrinum expression for DEGs that differed 
between the hybrid mef? genotypes by > 1.5-fold at one or 
more time points (Venables and Ripley 2010). NMDS is an 
ordination technique for dimension reduction that can be 
used on any measure of dissimilarity (i.e., distance). It uses 
an iterative regression-based approach to find a numerical 
solution for representing the distance relations among obser- 
vations from a high dimensional data set in k dimensions, 
where k is a small number relative to the dimensionality of 
the original data set. NMDS is nonmetric in the sense that it 
only attempts to maintain the rank order of distances in the 
spatial representation. The stress is a common statistic for 
assessing the fit of an NMDS solution (Venables and Ripley 
201 0) that can be interpreted as a measure of the distortion in 
the distance relationships that results from dimension reduc- 
tion. When conducting NMDS, we used 1 -ras a measure of 
dissimilarity where r is Pearson's correlation coefficient. This 
measure of dissimilarity ensures that all distances fall between 
zero and two and that the correlation across genes between 
samples determines the magnitude of dissimilarity, with pos- 
itively correlated samples treated as similar and negatively cor- 
related samples treated as dissimilar. We allowed a maximum 
of 150 iterations and conducted NMDS for k= 1-10 dimen- 
sions to examine how the stress changed as a function of the 
number of dimensions used for ordination. 

Results 

Effect of metl Genotype on Larval Growth and 
Metamorphic Timing in Backcrossed Hybrids 

Of the 30 individuals allowed to develop indefinitely, only two 



metf 



individuals failed to metamorphose within 365 



days and were thus scored as paedomorphs. The single 
recombinant for the mef7 region (ngfr ex/Att , rasdl mex/mex , 



metamorphosis at a relatively early time (1 64 DPH), consistent 
with metl associating more closely with map2k3 and ngfr 
than the other loci (discussed later). Twenty-seven individuals 
that were nonrecombinant for the metl region meta- 
morphosed. Among these, a highly statistically significant 
difference in metamorphic timing was observed between 



metf 



and metf' 



individuals (fig. 1/4). Despite this 



map2k3 n 



rail 



mex/mex 



and shmtr ex/mex ) completed 



pronounced difference in metamorphic timing, there was no 
difference between the mef7 genotypes in larval growth tra- 
jectory throughout the time period that tissues were collected 
(table 1 and fig. 16). These results show that mef 7 affects the 
timing of metamorphosis and the probability of expressing 
paedomorphosis without affecting early larval growth rate. 

Comparative and QTL Mapping of metl 

Previous studies identified ngfr as the closest protein-coding 
marker to mef 7 in the Ambystoma genetic linkage map (Smith 
et al. 2005; Voss and Smith 2005). To obtain linked markers 
for ngfr, we scanned genomic regions flanking ngfr in the 
Homo sapiens and Gallus gallus genomes and then selected 
loci that were good candidates for mef 7 . Th is strategy was not 
very successful because the ngfr chromosome in H. sapiens 
(HsalT) has undergone considerable rearrangement during 
evolution and the ngfr chromosome in G. gallus (Gga28) con- 
tains a ngfr-like gene sequence that does not show high se- 
quence identity to the Ambystoma ortholog. We eventually 
mapped map2k3 to within 5 cM of ngfr, an unexpected result 
because map2k3 is approximately 2.6 x 10 7 bases from ngfr 
in the human genome and these loci occur on different chro- 
mosomes in the chicken genome (fig. 2A). Linkage of 
map2k3-ngfr implicated a new set of candidate genes for 
mef7 and, in particular, genes from the Smith-Magenis syn- 
drome (SMS) genomic region (Hsa1 7p1 1 .2). Four additional 
loci (rasd7, rail , shmtl, and drg2) were subsequently mapped 
to establish linkage of SMS loci to map2k3. Several genes 
were then mapped within 5 cM on the opposite flank of 
ngfr, including sefd2, cllorfSl, ccm2, kif9, and aurkb. All 
of these genes except aurkb are syntenic within the same 
genomic scaffold in Xenopus tropicalis, and sefc/2 and ccm2 
are linked on chromosome 2 in G. gallus. However, the ngfr- 
setd2-cl lorf51-ccm2-kif9 conserved synteny group does not 
locate to the SMS conserved synteny group in the X. tropicalis 
genome, nor are these genes found together in fish or amni- 
ote genomes. These results show that ngfr and map2k3 flank 
the position of an ancestral chromosomal breakpoint that is 
not observed in representative anuran (Xenopus) and reptile 
(Gallus) genomes and therefore may be unique to Ambystoma 
and possibly other salamanders (fig. 24). 

The most likely position of mef7 was determined by QTL 
analysis using two different mapping panels: LAB (N= 90) and 
WILD2 (A/=496). In LAB, the segregation of discrete meta- 
morphic and paedomorphic phenotypes was observed to be 
largely consistent with a single gene model of inheritance 
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Fig. 1. — (A) Box plots comparing the ages at which spontaneous metamorphosis was completed by biphasic metf mex (n = 8) and metl 
(n= 19) hybrids. The difference in metamorphic timing between the alternate metl genotypes is highly statistically significant (one-sided Mann-Whitney 
test, 1/1/ = 142.5, P= 0.0002). (B) Scatter plot comparing the growth trajectories of metT ylex/mex hybrids (black circles) and met1 mex/Att hybrids (gray triangles). 
The six parameter dummy variable regression model described in Materials and Methods and table 1 was fit to the data. Stepwise deletion tests revealed that 
retaining a separate quadratic regression coefficient (F q qqs = 0.051, P=0.822), a separate linear regression coefficient (F qqqg = 0.858, P=0.356), and a 
separate intercept (F 112 o = 0.083, P= 0.774) for met1 mex/Att hybrids was not necessary. Thus, the curve shown corresponds to the selected model 
(F 2i i2i = 2,181, P« 0.0001) and implies that there is no difference between the growth trajectories of the alternate metl genotypes. 



Table 1 

Parameter Estimates, f Statistics, f Statistic P Values, and Adjusted R 2 Values for the Models Fit During the Selection Procedure Used to Compare 
the Growth Trajectories of r netT mx,mex and met1 mex/Att Hybrids 



Estimate 


Full Model 


Step 1 


Step 2 


Step 3 


B t ±SE 


-1.34±2.57x KT 1 


-1.30±1.96x 10"' 


-1.36±1.83x 10" 1 


— 1.36± 1.81 x 10" 1 


B t t value 


-5.19 


-6.63 


-7.46 


-7.52 


HB t ) = 0 


8.73 x 1(T 7 


1.03 x 10~ 9 


1.48 x 10" 11 


1.05 x 10" 11 


6 t ±SE 


-3.56 x 1(T 2 ±3.64x KT 1 


-1.11 x 10^1 1.47 x 10" 1 


1.57 x 10~ 2 ±5.47x 10~ 2 


Removed 


6 t t value 


-0.10 


-0.75 


0.29 


Removed 


HGt) = 0 


0.92 


0.45 


0.77 


Removed 


T,±SE 


1.21 x 10-'±8.15x 10~ 3 


1.20x 10^5.84 x10~ 3 


1.21 x 10~ 1 ±5.74x 10~ 3 


1.21 x 10~ 1 ±5.72x 10~ 3 


Tj t value 


14.88 


20.55 


21.09 


21.18 


HTi) = Q 


<2.00x 10~ 16 


<2.00x 10" 16 


<2.00x 10~ 16 


<2.00x 10" 16 


T 2 i ±SE 


-4.20 x 10~ 4 ±5.74x 10~ 5 


-4.11 x 10~ 4 ±4.06x 10~ s 


-4.12 x 10~ 4 ±4.05x 10" 5 


-4.12 x 10~ 4 ±4.04x 10~ 5 


T 2 , t value 


-7.32 


-10.13 


-10.16 


-10.21 


P(T 2 ,) = 0 


3.28 x 10" 11 


<2.00x 10~ 16 


<2.00x 10" 16 


<2.00x 10~ 16 


(GT) t ,±SE 


-7.44 x 10~ 4 ± 1.15 x 10~ 2 


1.81 x 10~ 3 ±1.96x 10~ 3 


Removed 


Removed 


(GT) t , t value 


-0.07 


0.93 


Removed 


Removed 


PKGT)t,] = 0 


0.95 


0.36 


Removed 


Removed 


(GT 2 ) t ,±SE 


1.83 x 10" 5 ±8.14x 10~ 5 


Removed 


Removed 


Removed 


(GT 2 ) t , t value 


0.23 


Removed 


Removed 


Removed 


P[(GT 2 ) t ,] = 0 


0.82 


Removed 


Removed 


Removed 


Adjusted ft 2 


0.9721 


0.9723 


0.9724 


0.9726 



Note. — To further emphasize the rationale used during the model selection process, additional decimal places are shown for adjusted R 2 relative to the other reported 
values. 



(Voss 1995). Accordingly, binary QTL mapping was used to 
map metl in LAB. A genome-wide scan for QTL using approx- 
imately 900 markers identified only one locus on linkage 
group 2 (LG2) with LOD scores exceeding an empirically 



determined significance threshold (LOD > 3.35, P=0.05; 
fig. 26). The maximum LOD score obtained by interval map- 
ping placed metl between conserved syntenies defined by 
ngfr-setd2-kif9 and map2k3-rasd1-rai1-shmt1 (fig. 2Q. 
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Fig. 2. — {A) Comparative map showing marker positions on Ambystoma linkage group 2 (AxLG2; center) and the positions of their orthologs in Gate 
gallus (left) and Xenopus tropicalis (right). (6) Results of a genome-wide scan for life history pathway QTL (metamorphosis vs. paedomorphosis) using the LAB 
mapping panel. A single significant LOD peak corresponding to met 7 was identified on Ambystoma LG2. (Q Local map of the met 7 genomic region on LG2 
based on the LAB mapping panel. (D) Result of a genome-wide scan for metamorphic timing QTL based on the WILD2 mapping panel. A single significant 
LOD peak corresponding to mef7 was identified on LG2. 



However, this position is based on a single recombinant ge- 
notype, and the LOD is equivalent when estimated at the po- 
sition of any of these flanking markers. To more accurately 
map metl, we turned to the larger WILD2 panel. In WILD2, 
we previously reported that -90% of individuals showed con- 
tinuous variation in metamorphic timing, and this variation 
was largely explained by metl (Voss and Smith 2005). Thus, 
we typed genes from the metl genomic region and per- 
formed QTL analysis of metamorphic timing to identify the 
most likely position for metl in WILD2. Gene orders were 
resolved to a finer degree using the WILD2 mapping panel 
(fig. 2D). The maximum LOD peak coincided with the position 
of an EST contig (nh_397) that is presumptively unique to A. 
tigrinum spp. because it does not show sequence identity to 
any NCBI nucleotide or protein-coding sequence. However, 
the LOD was approximately equivalent for a 3-cM interval 
defined by flanking genes aurkb and setd2/ngfr. The LOD 



score dropped approximately 2 LOD between ngfr and 
map2k3 and declined further for other SMS loci, although a 
minor peak was observed at the position of rasdl. Under the 
assumption that met 1 corresponds to a single gene, the most 
likely position for this gene is between aurkb and map2k3. 

Effect of metl Genotype on Gene Expression in 
Backcrossed Hybrids 

We identified 187 ZMMPSs and 6 EZMMPSs that measured 
differential transcript abundance between the alternate metl 
genotypes at one or more time points using limma (supple- 
mentary table S2, Supplementary Material online). Thus, there 
is no statistical evidence that this list of DEGs is enriched 
with EZMMPSs (odds ratio =1.41, lower 95% confidence 
limit = 0.47, upper 95% confidence limit = 3.59, P[odds 
ratio = 1] = 0.44). A few of these genes were identified as 
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Fig. 3. — Heat map showing the row scaled data from the 1 93 genes identified using limma. The columns correspond to individual GeneChips with the 
prefixes mm and mt corresponding to met1 mex/mex and met1 mex/Att hybrids, respectively. The time points in DPH are denoted by the numbers immediately 
following the lettered prefixes. Met1 mex/mex and met1 mex,Att GeneChips are separated by a vertical white line. The dendrogram to the left was obtained via 
hierarchical clustering of a Euclidean distance matrix. 



differentially expressed at the first five time points (42, 56, 70, 
84, and 98 DPH), but most were differentially expressed at the 
last time point (1 12 DPH). A substantial proportion of these 
DEGs exhibited temporal changes in both genotypes between 
84 and 98 DPH (fig. 3). However, between 98 and 1 12 DPH, 
most of these genes returned to expression levels similar to 
those observed between 42 and 84 DPH in met1 mex/mex indi- 
viduals (fig. 3). Conversely, the expression changes initiated 
between 84 and 98 DPH continued or intensified in met1 mex/ 
^individuals (fig. 3). The dummy variable quadratic regression 
approach identified 93 ZMMPSs and 3 EZMMPSs that 



measured different temporal patterns of expression between 
the alternate met! genotypes (supplementary table S3, 
Supplementary Material online). Hence, there is also no statis- 
tical support for the idea that the DEGs identified via maSigPro 
are statistically enriched with EZMMPSs (odds ratio =1.38, 
lower 95% confidence limit = 0.26, upper 95% confidence 
limit = 4.59, P[odds ratio= 1] = 0.49). Visual inspection of 
these 96 profiles identified 10 genes with robust expression 
differences between the alternate met! genotypes (table 2; 
supplementary fig. S1, Supplementary Material online). Two 
of the genes from this list map to LG2 (table 2), and 6 of the 1 3 
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Table 2 

The 10 Genes Identified via Dummy Variable Regression and 
Graphical Inspection with Temporally Consistent Differences in mRNA 
Abundance between Backcrossed Hybrids with Alternate met! 



Genotypes 


Gene 


Probe Set ID 


Sal-Site Version 3 ID 


Atp6v0c* 


SRV01135. 


_a 


_at 


Contig 18638 


Dhcr7 


SRV 00879 


_s_ 


at 


Contig40972 


Eif4a1 


SRV 00923. 


_a_ 


_at 


Mex_N M_001 41 6_Contig_9 


Gdi2 


SRV 00964 


_a_ 


_at 


Mex_NM_001 494_Contig_1 


H2afy2 


SRVJ3464 


a 


_at 


Tig_NM_018649_Contig_1 


Hdac2 


SRV 02393. 


a 


_at 


Contig40786 


Hsp90b1 


SRV01812. 


a 


at 


Mex_N M_003299_Contig_6 


NHO 


SRV_09160_ 


_a. 


_at 


Contig77873 


Raci* 


SRV 03102. 


_a. 


_at 


Contig20873 


Tomm70a 


SRV 03611. 


_a. 


_at 


Contig02244 



Note. — Sal-Site identifiers correspond to assembly version 3. Profiles of 
these genes are shown in supplementary figure S1, Supplementary Material 
online. NHO, no established human ortholog. Asterisk indicates genes that map 
to Ambystoma LG2. 



DEGs identified from both statistical approaches that map to 
LG2, located within 35 cM of metl. These patterns indicate 
that mef 7 genotype is associated with cis and trans transcrip- 
tional changes that commence late in the larval period, before 
morphological metamorphosis. 

Metl Is Associated with the Expression of Genes with 
Mitochondria] and Neurodevelopmental Functions 

Amphibian metamorphosis depends on transcriptional activa- 
tion of numerous biological processes that transform larval 
phenotypes into adult phenotypes (Shi 2000). Because our 
analyses suggested that the vast majority of metl's transcrip- 
tional effects are exerted late in the larval period, we used 
genes that were differentially expressed between met1 mex/ 
mex and met1 mex/Att hybrids at 1 12 DPH to identify biological 
processes associated with metl. As can be seen in table 3, 
DEGs upregulated in met1 mex/Att hybrids at 112 DPH were 
highly enriched for genes that function in the mitochondria, 
with particular enrichment for genes associated with electron 
transport (supplementary table S4, Supplementary Material 
online). The 24 DEGs that mapped to the mitochondrion 
ontology term are associated with a variety of functions in- 
cluding protein synthesis (mrpll3, mrpl30, and mrpsT), steroid 
metabolism {hint!), and apoptosis (eyes). In addition, six genes 
(ndufa4, ndufa7, ndufabl, ndufd, ucrc, and uqcrh) associ- 
ated with mitochondrial ATP synthesis coupled to electron 
transport were also upregulated in met1 mex,Att hybrids at 
112 DPH, indicating a change in mitochondrial energetics in 
individuals that metamorphosed (on average) at significantly 
earlier times. The DEGs that were upregulated in me ti mex/mex 
individuals (i.e., the genes downregulated in me t1 mex/Att indi- 
viduals) at 1 12 DPH were enriched with genes that mapped 
to a number of biological processes associated with the 
endomembrane system, translation, and neurodevelopment 



(table 3; supplementary table S5, Supplementary Material on- 
line). Of particular interest are several genes (actb, mtpn, rhoa, 
raci , rtn4, and ubc) associated with neurogenesis in mam- 
mals that were upregulated in met1 mex/mex hybrids relative 
to met1 mex/Att hybrids. Collectively, these results suggest 
divergent patterns of brain development and function be- 
tween m etl mex,mex and metl mex/Att hybrids. In particular, 
met jmex/Att hybrids increase transcription of genes associated 
with mitochondrial bioenergetics at a time that precedes 
metamorphosis. 



Expression of Hybrid DEGs in the Parent Species 

In a previous study, we characterized growth rates and 
neurodevelopmental gene expression in A. mexicanum and 
A. t. tigrinum larvae from 42 to 84 DPH (Page et al. 2010). 
In that study, A. t. tigrinum larvae grew at a faster rate, and all 
A. t. tigrinum showed definitive signs of morphological meta- 
morphosis by 84 DPH. Metamorphosis-associated genes were 
identified from A. t. tigrinum that changed in abundance as 
larval development proceeded, and these changes were par- 
ticularly pronounced at the later time points. In addition, 419 
probe sets with zero nucleotide mismatches between A. mex- 
icanum and A. t. tigrinum measured consistently different 
magnitudes of expression between species, with expression 
generally higher in A. t. tigrinum. Reexamination of how the 
193 DEGs identified from hybrids in this study using limma 
were expressed in the parent species revealed a robust pattern 
of gene expression divergence (fig. 4). The distribution of 
P values for the 772 r statistics we computed from the parental 
data differs greatly from the uniform distribution that would 
be expected if there were no expression differences between 
the species (fig. 5). This general pattern also holds for each 
time point when the P values are not pooled (not shown). 
Approximately 65% of these genes were differentially ex- 
pressed (a = 0.05) between A. mexicanum and A. t. tigrinum 
at one or both of the two earliest time points (42 or 56 DPH). 
Moreover, approximately 78% of these genes were differen- 
tially expressed between the parent species (a = 0.05) at one 
or more of the four time points that we sampled parental 
gene expression. These results show that many of the genes 
that are differentially expressed as a function of mef7 geno- 
type at 1 12 DPH in backcrossed hybrids are also differentially 
expressed between A. mexicanum and A. t. tigrinum, and that 
many of these differences are detectable in the parent species 
early during the larval period. 

To visualize the degree of similarity in expression between 
hybrid mef 7 genotypes and the parent species, we conducted 
NMDS on the 24 DEGs that were differentially expressed be- 
tween met1 m€ 



and metV 



hybrids by > 1.5-fold at 



one or more time points. The algorithm converged within 
150 iterations for solutions in /c=2-10 dimensions. As can 
be seen in figure 6A, the stress did not level off until k was 
beyond a number of dimensions that could be easily 
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Table 3 

Gene Ontology Terms that Were Significantly Enriched in Genes that Were Differentially Up and Downregulated in met1 mex,Att Hybrids Relative to 

met1 mexfim* Hy ^S 



GO ID 


GO Term 


N 


P 




Total 


Enrichment 


GO terms enriched 


in the DEGs upregulated in met1 mex,Att hybrids at 112 DPH 












0005743 


Mitochondrial inner membrane 


17 


2.32 x 10 


6 


58 


3.70 


0022900 


Electron transport chain 


9 


2.39 x 10 


5 


18 


6.24 


0005739 


Mitochondrion 


24 


2.17 x 10" 


4 


143 


2.12 


0006119 


Oxidative phosphorylation 


9 


2.72 x 10 


4 


24 


4.68 


\JKJ 1 JU / O 


M\/Hrr*Mpn ir*n trancmomhcano trancn^irtor ar+i\/i+v 
nyUJ UUCI 1 IUII UdllMll<_IIIUIcllltr Llall3fJUILd CIL.LI V 1 Ly 


3 


4 92 x 10" 


4 


21 


5.01 


UUUJ / J J 


^+ri i/-tr iral rr\nc+i+i lont r\T rihr*c/"imp 
-> LI ULIU 1 cl 1 L.UI IbLI LUL r 1 L <J 1 IIU'-OUIIIL. 


13 


9 44 x 10" 


4 


60 


2.85 


0042775 


Mi tnrhnnHrip 1 ATP c\/nthp*:i<: rounlpH plpctrnn tr;in<:nrirt 

1 VI 1 LULr IUI IUI lai 1 r jy 1 1 LI ICjIj LUU U ICLJ CICLII VJI 1 LI d 1 OLJL/I L 


5 


3 06 x 10 


3 


14 


5.35 


0016491 


Oxidoreductase activity 


13 


0.01 




79 


2.16 


001 5980 


Energy derivation by oxidation of organic compounds 


7 


0.02 




29 


3.01 


0006414 


Translational elongation 


10 


0.03 




57 


2.19 


0051234 


Establishment of localization 


25 


0.03 




213 


1.46 


0006120 


Mitochondrial electron transport, NADH to ubiquinone 


4 


0.04 




10 


4.99 


GO terms enriched 


in the DEGs downregulated in metr exlAn hybrids at 112 DPH 












0012505 


Endomembrane system 


14 


2.13 x 10 


4 


69 


3.10 


0008135 


Translation factor activity, nucleic acid binding 


6 


0.01 




22 


4.11 


0003924 


GTPase activity 


7 


0.01 




31 


3.41 


0022008 


Neurogenesis 


6 


0.02 




25 


3.59 


0006417 


Regulation of translation 


5 


0.03 




19 


3.94 



Note. — In total, 70 gene ontology terms were enriched in the list of DEGs that were upregulated in met1 msxJAtt hybrids and 32 gene ontology terms were enriched in 
the list of DEGs that were downregulated in metl msx/Att hybrids. The data in the table correspond to a subset of these results, which are shown in their entirety in 
supplementary tables S4 and S5, Supplementary Material online. P — modified Fisher's exact P value. N — the number of DEGs that annotate to a given ontology term. 
Total — the number of genes on the array that annotate to a given ontology term. Enrichment — the increase in DEGs annotating to a given ontology term relative to the 
number of genes expected given a random draw from the array of equal size. 



visualized. Nevertheless, there are several features of the 
NMDS solutions where k=2 (fig. 66) and k=3 (fig. 6Q 
that suggest there is: 1) hybrid intermediacy in gene expres- 
sion, 2) a change in both hybrid metl genotypes at 98 DPH 
that causes their expression profiles to more closely resemble 
A. t. tigrinum expression along at least one dimension, and 3) 
a reversal of the temporal trajectory initiated between 84 and 
98 DPH in m e ti mex/me x hybrids that is not observed in met1 mex/ 
Att hybrids (fig. 3). When k=3, these patterns are readily ap- 
parent in the values obtained for the horizontal axis depicted 
in figure 6C. The distributions of the coordinates for this di- 
mension are shown for several developmentally relevant 
groups in figure 6D. Collectively, our reexamination of 
hybrid DEGs in the parental species suggests that many 
differences between A. mexicanum and A. t. tigrinum in neu- 
rodevelopmental gene expression trace to genetic factors 
from the metl genomic region. 

Discussion 

Evolutionary History of the met 7 Genomic Region 

An important objective in evolutionary genetics is to identify 
genes for adaptive traits (Feder and Mitchell-Olds 2003). This 
objective has been pursued most often by testing a priori 



selected candidate genes (Voss et al. 2000, 2003), with 
fewer studies using unbiased approaches to identify candidate 
genes for QTL. A logical first step in searching for candidate 
genes is to cross-reference a QTL genomic region to a refer- 
ence genome or genetic map. In our case, some of the met1- 
associated genes that we initially mapped were not found to 
be syntenic in any vertebrate genome. Through trial and er- 
ror mapping of candidate loci and application of a functional 
genomics approach, we efficiently identified genes linked to 
metl . In the process, we obtained a sufficient number of loci 
to reconstruct the dynamic and lineage-specific history of 
chromosomal rearrangements that structured LG2 and 
the metl genomic region. LG2 primarily contains three ances- 
tral blocks of loci that correspond to ancestral vertebrate 
chromosomes (Voss et al. 2011). In the chicken genome, 
these blocks correspond to Gga12, Gga14, and a portion of 
chicken chromosome 2 (Gga2). Genes from these three 
blocks are not syntenic in Xenopus and thus were fused 
after the divergence of anurans and salamanders from a 
common ancestor. During salamander phylogenesis, intra- 
chromosomal rearrangements occurred between segments 
corresponding to Gga2 and Gga14, interleaving conserved 
syntenic groups of loci and creating unique gene orders. 
This is uncharacteristic of the Ambystoma genome, which 
has undergone a lower overall chromosomal rearrangement 
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Fig. 4. — Heat map showing row scaled data from the parent species for the 1 93 genes identified from backcrossed hybrids using limma. The columns 
correspond to individual GeneChips with the prefixes AM and AT corresponding to A. mexicanum and A. t. tigrinum samples, respectively. The time points in 
DPH are denoted by the numbers immediately following the underscore. GeneChips from different species are separated by a vertical white line. The 
dendrogram to the left was obtained via hierarchical clustering of a Euclidean distance matrix. 



rate than other vertebrates and especially mammals (Smith 
and Voss 2006). In fact, the chromosome associated with 
LG2 contains more lineage-specific, intrachromosomal rear- 
rangements than are predicted from any other Ambystoma 
linkage group (Voss et al. 201 1). Additional comparative map- 
ping studies are needed to determine whether this chromo- 
somal evolutionary history is unique to A tigrinum spp. or 
shared among all extant salamanders. Of particular relevance 
to this study is the fact that these chromosomal rearrange- 
ments brought several genes with associated functions into 
linkage during evolution. 



Gene Expression and the met 7 Genomic Region 

Our functional genomics approach accomplished two impor- 
tant objectives. First, by comparing the abundances of RNAs 
between individuals that inherited different metl genotypes, 
we efficiently identified expression differences for loci linked 
to metl. Second, by comparing transcription within develop- 
ing larval brains, we identified expression differences that are 
functionally associated with species-specific differences, in- 
cluding life history variation. Although it remains possible 
that metl corresponds to a single gene, the observation 
that several genes in linkage disequilibrium with meti were 
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Fig. 5. — Histogram showing the distribution of the 772 f statistic 
P values that were calculated from previously generated data on 
A. mexicanum and A. t. tigrinum. The dashed vertical line denotes 
P> |f| = 0.05. 



differentially expressed raises the possibility that more than 
one gene underlies the effect of metl on metamorphic 
timing and gene expression. Moreover, when considering 
genes mapped in this study, and other genes that are pre- 
dicted by synteny to map to the reconstructed metl region, 
we find that several sets of functionally associated genes were 
brought into linkage during evolution. These include genes 
associated with brain development (ngfr, ccm2, setd2, 
cspg5, and rasdl), transcriptional regulation (setd2, smarcd, 
spop, and mec/9), lipid biosynthesis (srebfl , scap, and pemt), 
and mitochondrial function (tmemll, atpafl, and nt5m). 
Notably, srebfl and scap encode proteins that physically inter- 
act to regulate transcription of genes associated with cellular 
lipid homeostasis (Brown and Goldstein 1999). Linkage of 
functionally associated genes suggests that metl presents 
characteristics of a supergene. A supergene is a group of 
linked loci that cosegregate alleles for specifying suites of 
coadapted traits (Joron et al. 2006). In the case of metl, 
genes that were brought into linkage may specify multiple 
brain phenotypic differences between metamorphic A. t. tigri- 
num and paedomorphic A mexicanum. It will be interesting to 
test the idea that metl is a supergene in natural populations 
where tiger salamanders express metamorphosis or paedo- 
morphosis as a polymorphism to exploit ephemeral and per- 
manent ponds, respectively. 

Although it is unknown which mechanism functions in 
A. mexicanum to yield a paedomorphic developmental out- 
come, a recent study showed that metl is a TH-responsive 
QTL (Voss et al. 2012). We extend this finding to hypothesize 
that metl affects the brain's response to TH during develop- 
ment, and as a result, this causes variation in transcription and 



timing of metamorphosis. We observed a large number of 
genes that were expressed differently between the hybrid 
metl genotypes by 112 DPH, and many of these DEGs 
were also differentially expressed between A. mexicanum 
and A. t. tigrinum as early as 42 DPH and extending to 84 
DPH. It is well established that TH activates new patterns of 
gene expression during anuran metamorphosis and vertebrate 
brain development (Shi 2000), including changes in cellular 
metabolism not unlike the highly enriched mitochondrial bio- 
logical processes that we observed in our study. In particular, 
genes that function in electron transport and ATP synthesis 
were enriched in our study, suggesting a critical role for oxi- 
dative phosphorylation in meeting the energetic demands as- 
sociated with neuronal differentiation and developmental 
processes that occur during metamorphosis. Overall, our re- 
sults show that high and low brain bioenergetics profiles are 
established very early in the larval period between A. t. tigri- 
num and A. mexicanum, but later in the case of hybrids. This 
may explain why these species show different larval growth 
rates and life history outcomes. 

Growth, Differentiation, and Metamorphic Timing 

The timing of metamorphosis is an ecologically and evolution- 
ary important trait because it is influenced by temporal as- 
pects of environmental factors that determine whether and 
when the adult habitat will become more suitable for growth 
and reproduction than the larval habitat. Consequently, the 
earliest and most influential ecological model of amphibian 
metamorphosis proposed that an organism's growth history 
is a critical determinant of metamorphic timing (Wilbur and 
Collins 1973). According to this model, there is a range of 
body sizes in which metamorphosis occurs, and within this 
range, metamorphosis will not proceed until the growth 
rate drops below a threshold value. Although there is strong 
empirical support for the idea that a threshold model is the 
appropriate framework for conceptualizing metamorphic 
timing (Morey and Reznick 2000), there is debate over 
whether growth rate or differentiation rate, two traits that 
are often highly correlated, is the crucial determinant 
(Smith-Gill and Berven 1979). At 112 DPH, we observed 
large-scale differences in transcription between metl geno- 
types that were not preceded by differences in growth rate. 
Moreover, many of these differences in gene expression 
appear to be associated with metamorphic differentiation 
events that are eminent in metl mexJAtt hybrids, but not 
metl mex/mex hybrids. Although the relevance of our findings 
remains to be established for natural populations, our results 
show that differences in metamorphic timing are not neces- 
sarily coupled to differences in larval growth history. 

Conclusion 

We showed in this study that metl locates to a region of the 
Ambystoma linkage map that has been structured by 
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EH H98 M112 Mex T112 Tig 

Dimension x 

Developmental group 

Fig. 6. — Results of the NMDS analysis: (A) scree plot showing the stress (in percent) for NMDS solutions based on different numbers of dimensions, (6) 
scatter plot of the coordinates obtained when k = 2, (Q scatter plot of the coordinates obtained when k = 3, and (D) box plots of various developmental 
groups' coordinates for the horizontal dimension shown in (O- In (6) and (0, green indicates A. mexicanum; blue, A. t. tigrinum; black, met1 mex/mex hybrids; 
red, met1 mex/Att hybrids. Circles = 42 DPH; triangles =56 DPH; + = 70 DPH; x =84 DPH; diamonds = 98 DPH; and inverted triangles = 1 12 DPH. In (D), 
EH = 42-84 DPH hybrids irrespective of genotype, H98 = 98 DPH hybrids irrespective of genotype, Mex = A mexicanum irrespective of sampling time, 
M 1 1 2 = 1 1 2 DPH met 7 mejl/mex hybrids, T1 1 2 = 1 1 2 DPH metl mex/Aa hybrids, and Tig = A. t. tigrinum irrespective of sampling time. 



chromosomal rearrangements, creating a gene order that 
has not been observed in any other vertebrates. Some of 
the loci in the met! region are associated with bioenerget- 
ics and brain development and function. Several of these 
genes, and many other genes that are not linked to metl, 
are differentially expressed between hybrids with alternate 
metl genotypes. In addition, we showed that many of the 
genes that were differentially expressed late in the larval 
period between hybrids with alternate metl genotypes 
were also differentially expressed between A mexicanum 
and A t, tigrinum at an earlier time of development. These 



results support the idea that differences between A mex- 
icanum and A t. tigrinum in metamorphic timing and neu- 
rodevelopmental gene expression are associated with 
genetic factors in the mef7 genomic region. Finally, we 
have shown that despite its profound effects on develop- 
mental timing and gene expression, met! does not influ- 
ence the larval growth trajectories of backcrossed hybrids 
prior to its effects on transcription in the brain. Overall, our 
study underscores the importance of using an integrative 
genomics approach to resolve genetic aspects of life history 
variation. 
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Supplementary Material 

Supplementary tables S1-S5 and figure S1 are available at 
Genome Biology and Evolution online (http://www.gbe. 
oxfordjournals.org/). 
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